Coarse-grained Interaction Potentials for Anisotropic Molecules 
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We have proposed an efficient parameterization method for a recent variant of the Gay-Berne 
potential for dissimilar and biaxial particles and demonstrated it for a set of small organic molecules. 
Compared to the previously proposed coarse-grained models, the new potential exhibits a superior 
performance in close contact and large distant interactions. The repercussions of thermal vibrations 
and elasticity has been studied through a statistical method. The study justifies that the potential 
of mean force is representable with the same functional form, extending the application of this 
coarse-grained description to a broader range of molecules. Moreover, the advantage of employing 
coarse-grained models over truncated atomistic summations with large distance cutoffs has been 
briefly studied. 
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I. INTRODUCTION 

The development of accurate, reliable and computa- 
tionally efficient interaction models is the main activity 
of molecular modeling. The need to attain larger simu- 
lated time scales and the excessive complexity of a wide 
range of molecular systems (e.g. biomolecular) has em- 
phasized the factor of computation efficiency as a domi- 
nant deliberation in choosing the appropriate interaction 
model for molecular simulations. In particular, grouping 
certain atoms into less detailed interaction sites, known 
as "coarse-graining", in one way of achieving such effi- 
ciency. 

Various coarse-grained (CG) approaches have been re- 
cently developed with such goal in mindi^*^. The im- 
plementation of coarse-graining models is usually divided 
into two distinct stages. The first is a partitioning of the 
system into the larger structural units while the second 
stage is the construction of an effective force field to de- 
scribe the interactions between the CG units. Typically, 
CG potentials of a pre-defined analytical form are param- 
eterized to produce average structural properties seen in 
atomistic simulations. Such analytical forms are chosen 
in a way to describe the governing interaction between 
the CG units^. The parameterizations are usually based 
on matching samples of potentials of mean forced, in- 
verse Monte Carlo data^ or certain atomistic potentials 
characteristicsS. The main concern of the present work 
is parameterizing a CG force field for the short-range 
attractive and repulsive interactions between ellipsoidal 
molecules and groups, based on atomistic potential sam- 
pling and potential of mean force. 

In molecular simulations, short-range attractive and 
repulsive interactions are typically represented using 
Lennard- Jones (6- 12) potentialsM: 



ULjir;i,j) = 4ej 



(1) 



where aij and are the effective heterogeneous inter- 
action radius and well-depth between particles of type 
i and j respectively and r is the inter-particle displace- 



ment. While the r~^ part has a physical origin in dis- 
persion or van der Waals interactions, the r^^^ repul- 
sion is chosen for mathematical convenience and is some- 
times replaced by exponential terms as well. For large 
molecules, the exact evaluation of the interaction po- 
tential of this type involves a computationally expensive 
double summation over the respective (atomic) interac- 
tion sites: 



U^ntiMl,M2) = V V Ua{nf,i,j) 



ieMi jeM2 



(2) 



where A4i and A^2 denote the interacting molecules and 
Ua{-) is the atomic interaction potential, e.g. Eq. QJ. In 
practice, a large distant interaction cutoff accompanied 
by a proper tapering is used to reduce the computation 
cost. More sophisticated and efficient summation meth- 
ods such as Ewald summation and the Method of Lights 
are also widely used^ . 

As an alternative approach. Gay and Bernei proposed 
a more complicated single-site CG interaction potential 
(in contrast to sophisticated summation techniques) for 
uniaxial rigid molecules which was generalized to dissim- 
ilar and biaxial particles later by Berardi et al as well^. 
We will refer to this potential as the biaxial-GB in the 
rest of this article. 

In response to the criticism of the unclear microscopic 
interpretation of the GB potentiaUS, we have recently 
used results from colloid scienceii to derive an approx- 
imate interaction potential based on the Hamaker the- 
oryi^ for mixtures of ellipsoids of arbitrary size and 
shape, namely the RE^ potential^^. Having a parame- 
ter space identical to that of Berardi, Fava and Zannoni-, 
the RE^ potential agrees significantly better with the nu- 
merically evaluated continuum approximation of Eq. |(2Jl, 
has no unphysical large distant limit and avoids the in- 
troduction of empirical adjustable parameters. 

In an anisotropic coarse-grained model, a molecule M. 
is treated like a rigid body. Neglecting the atomic details, 
each molecule is characterized by a center separation r 
and a transformation operator (a unitary matrix A or a 
unit quaternion q) describing its orientation. 
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In the first section of the article, we briefly introduce 
the RE^ potential followed by a review of the biaxial-GB 
potential and Buckingham(exp-6) atomistic model. The 
Buckingham(exp-6) potential is used in the MM3 force 
field^'^ and will serve as the atomistic model potential for 
parameterizations. The second section describes a pa- 
rameterization method which has been demonstrated for 
a few selected molecules, followed by an exemplar com- 
parison between the RE^ and the biaxial-GB potential. 
We will study the repercussion of internal vibrations, in 
contrast to the usually assumed proposition of the ideal 
stifFnessi*^, and propose an error analysis method to de- 
fine trust temperature regions for single site potentials. 
Finally, we will show that the potential of mean force 
(PMF) is representable with the same functional form 
for a wide range of temperatures. 



II. ATOMISTIC AND SINGLE-SITE 
NON-BONDED POTENTIALS 

In the MM3 force field^'', the van der Waals interac- 
tion is described in terms of Buckingham(exp-6) potential 
which is an exponential repulsive accompanied by a 
attractive term: 

C/™N;i, j) = e,, I^Ae-^--/'-- - C (^^^^ (3) 

where A, B and C are fixed empirical constants while (Ty 
and eij are heterogeneous interaction parameters specific 
to the interacting particles. Usually, Lorenz and Berth- 
elot averaging rules are used to define heterogeneous in- 
teraction parameters in terms of the homogeneous ones, 
i.e. atj = (ui -\- (Jj)/2 and tij = ^(-if-j- The hard core 
repulsion is usually described via a r~^^ term with an 
appropriate energy switching: 

vM^'\ny.^.2)^i[^f^ (4) 

where 7 is defined in a way to provide continuity at the 
switching distance. The interaction energy between two 
arbitrary molecules is trivially the pairwise double sum- 
mation over all of the interaction sites, i.e. Eq. (O. 

The dissimilar and biaxial Gay-Berne potential 
(biaxial-GB) is a widely used single-site model proposed 
by Berardi et al^ which is an extension of the original 
uniaxial description-'^ to biaxial molecules and heteroge- 
neous interactions. Based on the original Gay and Berne 
concept, the biaxial-GB is a shifted Lennard-Jones(6-12) 
interaction between two biaxial Gaussian distribution 
of interacting sites. In this coarse-grained model, each 
molecule is described by two diagonal characteristic 
tensors (in the principal basis of the molecule) S and E, 
representing the half radii of the molecule and the 
strength of the pole contact interactions, respectively. 
As mentioned earlier, the orientation of a molecule is 



described by a center separation vector r and a unitary 
operator A, revolving the lab frame to the principal 
frame of the molecule. 

The biaxial-GB description for the interaction between 
two molecules with a center separation of V\2 = r2 ^ ri 
and respective orientation tensors Ai and A2 is defined 
as: 

C/f,^(ri2,Ai,A2) = 



/ \ 12 




{ ^ 




\hi2 +(JcJ 


\hi2 +(JcJ 



where Cq and are the energy and length scales, 7712 and 
X12 are purely orientation dependant terms^ and hi2 is 
the the least contact distance between the two ellipsoids 
which are defined by the diagonal covariance tensor of 
the assumed Gaussian distributions. The orientation 
dependant terms (7712 and X12) describe the anisotropy 
of the molecules. 

We have recently proposed a single-site potential, 
namely RE^ — giving the approximate interaction energy 
between two hard ellipsoids in contrast to the tradition of 
the Gaussian clouds, initiated by Gay and Berne^. The 
orientation dependence of the RE^ potential fall at large 
distances, reducing asymptotically to the interaction en- 
ergy of two spheres. Moreover, it gives a more realistic 
intermediate and close contact interaction using a heuris- 
tic interpolation of the Deryaguin expansiorki^si^. Being 
a shifted Lennard-Jones(6-12) potential, the biaxial-GB 
fails to exhibit the correct functional behavior for large 
molecules i'^. The attractive and repulsive contributions 
of the RE^ potential are respectively: 

(Ai, A2, ri2) = (1 + 3^712X12^) X 

n n [^Stt^ ^'^^ 

(Ai,A2,r,2) = ^(^) [^+Y,m2Xi.j^J^ 

n n ( (.)/f\poO (6b) 

where Ayi is the Hamaker constant (the energy scale), 
CTc is the atomic interaction radius and a^' , a^^ and ctI*'' 
are the half-radii of ith ellipsoid (i=l,2). The terms 7712, 
X12 and h\2 are defined in parallel to the biaxial-GB 
model and thus, are described in terms of the same 
characteristic tensors. 

The structure tensor Si and the relative potential well 
depth tensor E; are diagonal in the principal basis of zth 
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molecule and are defined as: 



S,=diag{a«,aW,a«} 



(7a) 



where si*-* , Ey'' and i?!'' are dimensionless energy scales 
inversely proportional to the potential well depths of the 
respective orthogonal configurations of the interacting 
molecules (aa, bb and cc, Table P). For large molecules 
with uniform constructions, it has been showr>i^ that the 
energy parameteres are approximately representable in 
terms of the local contact curvatures using the Deryaguin 
expansioni^: 



diag{£;«,i?W,£;W} (7b) 



The most general parameter space of the RE^ potential 
contains the energy and length scales, the characteristic 
tensors and the parameters specifying the relative ori- 
entation of the ellipsoids to the molecules. In order to 
overcome the degeneracy of the parameter space and to 
guarantee the rapid convergence of the optimization rou- 
tines, a two-stage parameterization is proposed. In the 
first stage, the center and principal frame of the molecule 
will be fixed at the centroid and the eigenbasis of the in- 
ertia tensor. A preliminary optimization in the reduced 
parameter space yields to an approximate parameteriza- 
tion. In the second stage, the results of the first stage will 
be taken as the initial guess, followed by an optimization 
in the unconstrained variable space. This two-stage pa- 
rameterization will theoretically result in superior results 
for molecules with imperfect symmetries. 



Ei = CTcdiag 



CTuCTr ax(Jz crx(Tv 



(8) 



B. Sampling and Optimization 



The assumptions leading to these estimations are not 
valid for the studied small organic molecules. There- 
fore, we will cease to impose further suppositions and 
take these three scales as independent characteristics of 
a biaxial molecule. Computable expressions for the ori- 
entation dependent factors of the RE^ potential (ryi2 and 
X12) among with the Gay-Berne approximation for hi2 
has been given in the Appendix 



III. PARAMETERIZATION FOR ARBITRARY 
MOLECULES 



The Principal Basis and The Effective Center of 
Interaction 



Associating a biaxial ellipsoid to an arbitrary molecule, 
one must define an appropriate principal basis and a 
center of interaction for it beforehand, according to the 
used coarse-grained model. Although there's no trivial 
solution to this problem, the centroid and the eigenba- 
sis of the geometrical inertia tensor of the molecule are 
promising candidates and may be taken as suitable ini- 
tial guesses as they yield to the correct solution at least 
for the molecules with perfect symmetry. For a molecule 
consisting of N particles, the centroid is defined as: 



Vr = 



TV 



(9) 



and the principal basis is the eigenbasis of the geometrical 
inertia tensor Ig given by: 



JV 



l9 = II(^?l-^i'^^*) 



(10) 



where is the position of ith atom. 



Physical and symmetrical considerations lead to the 
proposition that a sampling of the pole contact interac- 
tions between two biaxial particles is essentially sufficient 
to reproduce the interaction for all configurations. There 
are 18 different orthogonal approaching configurations 
(pole contacts) between two dissimilar and biaxial par- 
ticles (Table nj. Based on physical grounds, we optimize 
the parameter space for the important characteristics of 
the sampled orthogonal energy profiles, i.e. potential 
well depth, potential well distance, the width of well 
at half depth and the soft contact distance. This 
parameterization fashion is guaranteed to produce a 
satisfactory reconstruction of the most crucial region of 
interaction. 

The geometry of the molecules were initially optimized 
using TINKER molecular modeling package^^^ with the 
MM3 force field. We have used the same force field 
to sample the interaction energy for the orthogonal 
configurations. 

Given a parameter tuple p, we denote the potential 
well depth, well distance, well width at half depth and 
the soft contact distance for ith orthogonal configuration 
predicted by the RE^ potential by ?7m(«;p), i?m(«;p), 
iy(i;p) and Rsc{i',p) respectively. The same potential 
well specifications calculated from the atomistic sum is 
denoted by scripted letters. An appropriate cost function 



1 "' 

o(p)^^i: 



wr„ 



i=l 

R7n{i;p) - Tlm{i] 



WW 



W{i;p)-W{i) 
Wo 



n 



(11) 
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where fJo is a normalization factor: 

no - +^-^ + <J^ E^"'''"^'^- (12) 

We have chosen Uq as min and TZq, Wq and 

7?.scO as min {yV„i(i)} based on physical considerations. 

is the number of orthogonal profiles (12 and 18 for ho- 
mogeneous and heterogeneous interactions respectively) 
and (uif/, Wfl, wvy, Wii;,^) are fixed error partitioning 
factors for different terms, set to (1.0,3.0,2.5,1.0) 
in order to emphasize on the structural details. We 
have also included a fixed error weighting according 
to the Boltzmann probability of the appearance of the 
corresponding profiles. One expects higher amplitude of 
relative appearance for orientations with deeper wells, 
which justifies the requisite of higher contribution in the 
cost function. We have also chosen jS as l/(|Wm(i)|) in 
order to avoid deep submergence of the lower energy 
orientations. 

Further implications such as matching the large 
distance behavior will be regarded as constraints on 
the parameter space, leaving the defined cost function 
unchanged. 

The nonlinear optimization procedure consists of a 
preliminary Nelder-Mead Simplex search followed by 
a quasi-Newton search with BFGS Hessian updatesii. 
The whole parameterization routine is coded in MAT- 
LAB/Octave and is freely available^®. The procedures 
of sampling and parameterization are purely automated 
and requires only a Cartesian input file. The interac- 
tion parameters for the homogeneous interaction of a set 
of small prolate and oblate organic molecules has been 
provided in Table iQlJ. It is noticed that the provided 
half radii agree significantly better with the molecular 
dimensions compared to the biaxial-GB parameteriza- 
tiona^, reflecting the precise microscopic interpretation 
of the RE^ potential. 



C. Large Distance Analysis 

The cost function defined in the previous section fo- 
cuses on close contact regions only. In order to achieve 
the correct large distant limit as well, we will constrain 
the variable space by matching the asymptotic behavior 
of the RE^ potential with the atomistic summation. The 
asymptotic behavior of the RE^ potential is described as: 

16 

lim r°2C/flB2(ri2,Ai,A2) = -— A12 det[Si] det[S2] 

ri2— foo y 

(13) 

The atomistic summation defined by Eq. (|2Jl and Eq. Q 
exhibits the same asymptotic behavior, which together 



with Eq. (|13|l results in the relation: 

A12 det[Si] det[S2] - ^ E E ^^^4 (14) 

The summation appearing in right hand side is most eas- 
ily evaluated by a direct force field parameter lookup. 
Applying such a constraint guarantees the expected large 
distant behavior while leads to a faster parameterization, 
reducing the dimensions of the variable space. Uncon- 
strained optimization routines are still applicable as one 
may solve Eq. (|14|l for A12 explicitly. A graphical com- 
parison between the biaxial-GB and the RE^ potential 
has been given for the homogeneous interaction of the 
pair PeryleneiSi has been sketched in Fig. ^ The large 
distance convergence of the RE^ potential is noticed in 
contrast to the divergent behavior of the biaxial-GB po- 
tential, which is due to the non-vanishing orientation de- 
pendent pre-factors. Although the energy contribution is 
small at this limit, it is not generally negligible, e.g. the 
large distant separability of the orientation dependence 
of the model potential alters the nature of the phase di- 
agram and the long range order of a hard rod fluid in 
general2£. 



D. Heterogeneous Interactions 

The heterogeneous interaction between two molecules 
Ml and M2 is calculable by equations and (|6b|) 
once the characteristic tensors of each molecule (S and E) 
along with the heterogeneous Hamaker constant AmiM2 
and the atomic potential radius <tmiM2 are available. 
The heterogeneous Hamaker constant may be evaluated 
directly with a force-field parameter lookup. Moreover, 
the arithmetic mean of (tmiMi and aM2M2 is a rea- 
sonable estimate for the heterogeneous interaction ra- 
dius, crMiM2- Therefore, the homogeneous interaction 
parameters of the molecules Aii and A^2 are sufficient 
to describe their respective heterogeneous interaction us- 
ing the RE^ potential. Apparently, there is no trivial 
mixing rule available for the energy scale of the biaxial- 
GB potential. Inspired by the atomic mixing rules and 
the theory of the Gay-Berne potential, we have used 
Berthelot's geometric averaging rule for this purpose. An 
instance of a heterogeneous interaction has been illus- 
trated in Fig. ||2Jl for the pair Perylene (oblate) and Sex- 
ithiophene (prolate)^^. The results are quite promising 
for a coarse-grained model; However, further optimiza- 
tion will theoretically yield to superior results. Conclud- 
ing from the graphs, the RE^ potential performs signif- 
icantly better at end-to-end and cross interactions com- 
pared to the biaxial-GB. The error measures (ilfl£;2 — 
6.5 X 10"'^, r^GB = 7.7 X 10"'^) agree with this observa- 
tion. 
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E. The advantages over practical atomistic 
implementations 

As mentioned before, the atomistic evaluation of long- 
range interaction potentials involve computationally ex- 
pensive double summations over the interaction sites, re- 
sulting in a quadratic time cost with respect to the av- 
erage number of interactions sites. However, the average 
computation time of a single-site potential is intrinsically 
constant, regardless of the number of interacting atoms. 
These observations have been quantified in Fig. (^J which 
is a comparison between the computation time of an ex- 
act LJ(6-12) atomistic summation and an efficient imple- 
mentation of the RE^ potential^* . Concluding from the 
graph, employing the RE^ potential for molecules con- 
sisting as low as ~ 5 atoms (or ~ 25 overall atomic in- 
teractions) is economic. 

The atomistic summations are practically employed 
with a proper large distance atomic cutoff in order to 
reduce the computation time. In the presence of large 
distance cutoffs, long range correction potential terms^ 
are usually used to compensate the submergence of par- 
ticles beyond the cutoff distance. 

Considerable errors may be introduced by choosing 
small atomic cutoff distances compared to the dimensions 
of the interacting molecules. Therefore, it is expectable 
that a CG model yield to relatively better results com- 
pared to truncated atomistic summations in certain con- 
figurations. A figurative situation is the end-to-end inter- 
action of two long prolate molecules. In such configura- 
tions, usual atomic cutoffs (~ 2.5a) can be small enough 
to dismiss the interaction between the far ends of the 
molecules. Moreover, long range correction terms are of 
little application in this case due to the excessive inho- 
mogeneity and the small number of interacting particles. 

This effect has been illustrated in Fig. for Pen- 
tacene moleculei^. The first panel is a semi-log plot of the 
relative error for the RE^ potential together with three 
atomistic approximations with different cutoffs (6, 9 and 
12 A). The discontinuity of the truncated atomistic sum- 
mations is a result of hard cutoffs. In a proper atomic im- 
plementation, tapering functions are used to avoid such 
discontinuities. Concluding from the graphs, the CG de- 
scription introduces less error in all ranges of this configu- 
ration compared to the truncated atomistic summations, 
even with unusually large atomic cutoff distance (12 A). 
Moreover, the evaluation of the CG interaction potential 
requires a considerably lower computation time. 



IV. INTERMOLECULAR VIBRATIONS AND 
SINGLE-SITE POTENTIALS 

In this section, we study the proposition of ideal rigid- 
ity of the molecules, which is widely assumed in single- 
site approximations of extended molecules, including our 
own study in the previous sections. The samplings are 
usually taken from the relative orientations of the un- 



perturbed and geometrically optimized structures. The 
resulting parameterization will be used in molecular dy- 
namics simulations in which internal vibrations may not 
be negligible. We will introduce a method to estimate the 
error introduced by this supposition in the first part of 
this study. A parameterization based on the Potential of 
Mean Force (PMF) is probably the best one can achieve 
with the coarse-grained models, although the samplings 
are expensive. We will study such parameterizations in 
the second part. 

A. Analysis of the Mean Relative Error 

The PMF for the interaction of semi-rigid molecules 
in an arbitrary ensemble may be expressed as an addi- 
tive correction term to the the interaction potential of 
the respective rigid molecules. We will show that these 
correction terms are expressible in terms of statistical 
geometric properties of a molecule in the ensemble. The 
PMF between two molecules A4i and A^2 with a mean 
center separation of ri2 = r2 — ri and mean orientation 
tensors Ai and A2 is defined as: 

[/p„/(ri2,Ai,A2) - (C/(Xi,M2)) (15) 

where (•) denotes the ensemble averaging. The mean lo- 
cation of intermolecular particles are expected to remain 
unchanged compared to the unperturbed structures for 
a large range of temperatures as the internal structures 
of semi-rigid molecules are mainly governed by harmonic 
bond stretching and angle bending potentials. 

We may assume the location of each particle as a ran- 
dom variable, sharply peaked at its mean value. There- 
fore, we denote the location of ith particle measured in 
its principal coordinate system by: 

Ti = f; + SVi (16) 

where Fj = (rj) and 6ri is a displacement due to internal 
vibrations with vanishing average. The PMF of the in- 
teraction between the molecules Aii and A^2 is defined 
as: 

C/p„/(ri2, Ai, A2) = / ^'idl'^i-'^J II 

(17) 

where Ua{-) is the atomistic interaction potential. It is 
easy to show that up to the second moments: 

(||r,-r,.||)~||r».|| + 

^i:Va,(.,-..,,..,(l-<|^) ,18, 

where r^^ = f ^ — ij . We have neglected the covariance be- 
tween the coordinates. Using the last relation, we reach 
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to a second-order estimate of Eq. (|17|l : 

f/pm/(ri2, Ai, A2; A^l, Al2) ^ 

J2 E f/a(||f.-r,||;*,j) + 

ieAii jeAi2 \ fe=i 

(4.e,)^C/^(||rO||;^,J) ^ 



^0 112 



EVar(<5r. - ^r,).efc || ; (19) 

fc=i " / 

The first term of the right hand side is the interaction 
energy of the averaged structures, where the remaining 
terms are second-order corrections. The RMS of the rel- 
ative error introduced by neglecting the correction terms 
(e.g. the error in parameterizations based on unper- 
turbed samplings) may be evaluated formally via the fol- 
lowing integral: 



1 



exp 



5U{uj) 



(20) 

where w is a relative orientation, dN^ is a differential 
measure of orientations near w, ^{T) being the ensem- 
ble and 5U [lo) is the second-order correction defined by 
Eq. p9|) . £o{T) is the normalization factor defined as: 



exp 



en(T) 



U{u) 
ksT 



dK: 



(21) 



In practice, the spatial variance of each particle 
in an ensemble is most easily obtainable through an 
MD simulation. Once the statistical information are 
accessible, £{T) is most easily evaluated by Monte Carlo 
integration. Neglecting the covariance between the 
dislocation of the particles, we are implicitly overlooking 
the stretching and bending of the molecules at close 
contact configurations. Although our proposed error 
analysis disregards this phenomenon, it still measures 
the introduced error due to purely thermal vibrations. 

We have demonstrated this error analysis method for 
three different molecules in a large range of temperatures. 
The statistical information was extracted from several 
MD simulation snapshots with the aid of TINKER molec- 
ular modeling packagaiSi, each with 32 molecules and 
with periodic boundary conditions in an NVT ensemble 
(Fig. [SJ. For each isothermal ensemble, the RMS error 
has been evaluated using the MC integration of Eq. H2U|) 
for 10^ random orientations. The relation between the 
RMS error and the temperature is noticeably linear. The 
linear regression analysis has been given at Table (jIII|l . 
According to the required degree of precision, one can de- 
fine a trust region for the temperature using diagrams like 



Fig. (0). For example, a mean relative error of less than 
10% is expected for temperatures less than 1500K in a 
homogeneous ensemble of Benzene molecules, concluding 
from the graph. It is also concluded that the studied pro- 
late molecule (Sexithiophene) exhibits a higher relative 
error due to its considerably higher elasticity, compared 
to the oblate molecules (Perylene and Benzene). 



B. Parameterizations based on the Potential of 
Mean Force 

The error introduced by the assumption of ideal 
rigidity may not be negligible for certain purposes, 
concluding from the previous analysis. However, a 
coarse-grained potential which is parameterized on a 
PMF basis is theoretically advantageous as it is expected 
to describe the mean behaviors closer to the atomistic 
model. Phenomenologically speaking, the internal 
degrees of freedom will soften the repulsions at close 
contacts while the thermal vibrations are expected to 
smoothen the orientation and separation dependencies 
of the interaction. 

The PMF for a given macroscopic orientation may be 
evaluated through a Constrained Molecular Dynamics 
simulation (CMD) process with appropriate restrains. 
We have used harmonic restraining potentials for the 
center separation vector and on the deviations from the 
desirable principal basis for each molecule in order to 
keep them at the desired orientation. In order to reduce 
the random noise of the evaluated PMF, we applied a 
fifth-order Savitsky-Golay smoothing filter followed by a 
piecewise cubic Hermite interpolating polynomial fitting 
to the PMF samples. 

Fig. © is a plot of the evaluated PMF between the 
pair perylene for the cross configuration be. The upper 
(and interior) plots refer to higher temperatures. The 
expansion of the potential well width at lower tempera- 
tures is related to the tendency of the molecules to bend 
and stretch and thus, resulting in a softer interaction 
while the shift of the soft contact and potential well dis- 
tance along with the elevation of the potential well is 
associated to the thermal vibrations and hence, the ex- 
pansion of the effective volume of the molecules. The 
temperature-dependant parameterizations (based on the 
evaluated PMF) given at Table (|1V() justifies these qual- 
itative discussions. One may associate the contraction of 
the molecule {a^, <Jy and <Tz) and the expansion of the 
atomic interaction radius at lower temperatures to con- 
traction of the molecule at rough repulsions and widen- 
ing of the potential well, respectively. Furthermore, the 
expansion of the molecule volume at higher tempera- 
tures reflect the overcoming of thermal vibrations to the 
flexibility of the molecule. According to tables (QJ and 
(|IV|) . the overall error measures (il) for PMF parameter- 
izations closely match the same measure for the unper- 
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turbed structure. Thus, the same functional form may 
be used to represent the PMF as well. 



where B12 is defined in terms of the orientation tensors 
Ai and relative well-depth tensors E^: 



V. CONCLUSION 

We have proposed and demonstrated a parameter- 
ization method for the RE^ anisotropic single-site in- 
teraction potential which leads to a globally valid de- 
scription of the attractive and repulsive interaction be- 
tween arbitrary molecules. Unlike the biaxial-GB'^, the 
RE^ potential gives the correct large distant interac- 
tion (Fig. ^ while having a superior performance in the 
close contact region (Fig. The Potential of Mean 
Force is representable with the same functional form of 
the RE^ potential. Compared to the parameterizations 
given a.iim, The structure tensors agree significantly better 
with the spatial distribution of the intermolecular parti- 
cles. It has also been shown that the coarse-gained mod- 
els perform significantly better at certain configurations 
in comparison to a truncated atomistic summation with 
large distance cutoffs. 
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APPENDIX A: THE ORIENTATION 
DEPENDENT TERMS 

We will briefly quote computable expressions for the 
orientation dependant terms from the original articled. 
The term xi2 quantifies the strength of interaction with 
respect to the local atomic interaction strength of the 
molecules and is defined as: 



Bi2(Ai,Ai) = AfEiAi+A^EsAa. (A2) 

The term 7712 describes the effect of contact curvatures 
of the molecules in the strength of the interaction and is 
defined as: 



r?i2(Ai, A2,ri2) 



det[Si]/(T? + det[S2]/fT: 



[det[Hi2]/(ai+CT2)] 



TtI' (A3) 



The projected radius of zth ellipsoid along fi2 (di) and 
the tensor H12 are defined respectively as: 



f7,(A,,fi2) = (rf2A/ ^Aifi2 



-1/2 



(A4) 



Hi2(Ai, A2,fi2) = — AfS?Ai + — Ai^S^A2. (A5) 

(Jl (72 



There's no general solution to the least contact distance 
between two arbitrary ellipsoids (/112). The Gay-Berne 
approximation^'^" is usually employed due to its low com- 
plexity and promising performance: 

h?2'' = \\ri2\\ - (A6) 
The anisotropic distance function tTi2^ is defined as: 



ai2 = ( hj,G-,Hu) ' (A7) 



where the symmetric overlap tensor G12 is: 



G12 = A^SjAi + AjS^Aa (A8) 



Xl2(Ai, A2,fl2) = 2ff2Bi2^(Ai, A2)fl2 



(Al) 
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squared. 

The chemical fornmla for the studied molecules are: 
(1) Perylene [C20H12] (2) Pyrene [CieHio] (3) Coroncne 
[C24H12] (4) Benzene [CeHs] (5) Sexithiophene [S6C24H14] 
(6) Pentacene [C20H16] (7) Anthracene [C14H10] (8) Naph- 
thalene [CioHs] (9) Toluene [CrHg]. 
J. A. Cuesta and D. Prenkel, Phys. Rev. A 42 (1990). 
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TABLE I: The 18 orthogonal configurations of two dissimilar biaxial particles. A unitary operator (U) followed by a translation 
is applied to the second particle to reach the desired configuration. We adopt the naming scheme introduced by Berardi et aii. 
The operator Re denotes a 7r/2 rotation with respect to the axis e. A two-letter code is attached to each configuration with 
respect to the faces perpendicular to connecting vector of the ellipsoids. A prime is added if one or three axes are antiparallel. 
Italic codes refer to configurations which are degenerate in homogeneous interactions. 



u 


ri2 ea;i 


ri2 e^i 


ri2 


I 


aa 


bb 


CO 




ab 


ha' 


cc' 


Ry 


ac' 


bb' 


ca 


R-a; R-y 


ac 


ba 


cb 


R-z R-x -f^y 


aa' 


be' 


cb' 


R^R^ 


ab 


be 


ca 



TABLE IL RE^ potential parameters for homogeneous interactions of selected molecules—. The oblate molecules are: (1) 
Perylene (2) Pyrene (3) Coronene (4) Benzene. The prolate molecules are: (5) Sexithiophene (6) Pentacene (7) Anthracene (8) 
Naphthalene (9) Toluene. 



Mol. No. 


yli2(10^ Kcal/mol) 


<7c(A) 


a. (A) 


CJy{k) 


a. (A) 


Ex 


Ey 


E, 


0(10-=^) 


Oblate: 




















(1) 


36.36 


3.90 


4.20 


3.12 


0.49 


3.96 


2.39 


0.49 


9.6 


(2) 


28.36 


3.91 


4.24 


3.07 


0.45 


3.98 


2.35 


0.43 


9.9 


(3) 


21.01 


3.83 


4.27 


4.26 


0.54 


2.84 


2.85 


0.35 


12.5 


(4) 


84.95 


3.99 


2.14 


1.82 


0.36 


4.60 


3.70 


1.03 


6.8 


Prolate: 




















(5) 


49.44 


4.07 


10.96 


1.99 


0.46 


6.30 


1.16 


0.35 


10.7 


(6) 


37.46 


3.88 


6.56 


2.28 


0.47 


5.52 


1.65 


0.43 


9.7 


(7) 


45.51 


3.85 


4.19 


2.25 


0.44 


6.83 


2.48 


0.62 


9.2 


(8) 


37.76 


3.82 


3.09 


2.20 


0.49 


4.59 


3.39 


0.77 


10.9 


(9) 


23.19 


3.75 


2.72 


2.04 


0.57 


4.61 


3.29 


1.00 


10.6 



TABLE IIP Linear regression analysis between {SU/U)rms and T for a few selected moleculesifii. The molecules are: (1) 
Benzene (2) Perylene (3) Sexithiophene. The linear relationship is defined as {5U/U)rms ~ AT+B in all cases. 

Mol. No. AjlQ-^K'^) B(10~^) 



(1) 0.66 1.0 0.987 

(2) 1.39 -4.1 0.983 

(3) 2.35 1.3 0.995 



TABLE IV: RE'^ potential parameters for the homogeneous interactions of the pair perylene at different temperatures. 



Temperature (K) 


yli2(10'' Kcal/mol) 


ae(A) 


(Tx{k) 


ay (A) 




Ex 


Ey 


E, 




100 


14.59 


4.83 


3.77 


2.69 


0.10 


3.56 


2.13 


0.41 


15.2 


300 


2.04 


4.51 


3.97 


2.85 


0.24 


3.52 


2.18 


0.43 


14.1 


500 


1.21 


4.40 


4.04 


2.91 


0.30 


3.49 


2.20 


0.44 


12.6 


700 


0.99 


4.33 


4.09 


2.97 


0.32 


3.53 


2.19 


0.44 


12.4 


900 


0.82 


4.31 


4.09 


2.98 


0.35 


3.43 


2.15 


0.44 


12.6 
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FIG. 1: A comparison between the RE and biaxial-GB potentials for the homogeneous interaction of the pair perylene in a 
set of uniform random center separations in the range [SA, 50A] along with random rotations. The Gay-Berne approximation 
has been used for the least contact distance, (a) A log-log plot of Uhe^ /Umm3 against Umms (Mean=-0.002, SD=0.08) (b) A 
log-log plot of Ugb/Umm3 against Umms (Mean=-0.87, SD=0.54) 
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FIG. 2: The heterogeneous interaction between the pair perylene (oblate) and sexithiophene (prolate) for the 18 orthogonal 
configurations. The black thick lines denote the RE^ potential, the red dashed lines refer to the biaxial-GE^ and the reference 
atomistic summation (MM3^'') is denoted by blue thin lines. A combination of homogeneous interaction parameters f Table HB 
have been used without further optimization. The error measures are: il^E^ = 6.5 x 10~^, Qgb = 7.7 x 10~^. The graphs are 
grouped in five plates as: side-by-side (A), cross (B), T-shaped 1 (C), T-shaped 2 (D) and end-to-end (E) interactions and are 
labeled according to the notation introduced in Table (0. 



12 



E 2.5 



.2 2 

i. 1.5 

E 
o 

O 1 



Interaction Sites 



FIG. 3: The average computation time of exact LJ(6-12) atomistic summation with respect to the average number of interacting 
sites (blue dashed line) and RE^ single-site potential (red continuous line). 




FIG. 4: A comparison between the truncated atomistic descriptions with hard atomic cutoffs and the RE'^ potential for the end- 
to-end interaction of the pair Pentacene. Thick lines denote the RE^ potential while thin lines represent atomistic summations 
with different atomic cutoffs (6, 9 and 12 A). (A) Logarithmic relative error of the RE^ potential and truncated atomistic 
summations (with respect to the exact atomistic summations) vs. center separation. (B) Time consumption of different 
approximations vs. center separation. 
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FIG. 5: Relative deviations from the PMF vs. temperature for three different molecules. The signs indicate the MD simulation 
data. The continuous lines are linear regressions. (1) Plus signs: Benzene (2) Cross signs: Perylene (3) Dots: Sexithiophene. 
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FIG. 6: Potential of Mean Force between the pair perylene for the cross configuration be. The dashed line indicate the 
interaction potential of the unperturbed structures while the continuous lines, ordered descending with respect to their well- 
depths, represent the PMF at temperatures lOOK, 300K, 500K, 700K and 900K, respectively. 



